% Scripts to calculate the LT dependency of prompt penetration
% By correlating IEF with EEF from Jicamarca ISR
% While similar scripts exist in jic_julia_delh_proc_2001_20XXXXX file,
% this is the updated one.
%% LT dependency
%
clear;
if ~exist('ace_all','var'),
    
    if ispc == 1,
        
        load E:\projects\ace_tensor\acedata\ace_2000_2010.mat ace_all;
        ace_all(:,1) = ace_all(:,1)+(17)/(60*24);%Try to advance 60 VS model/ ace time minutes
        load('C:\Manoj\projects\ace\jcamarca_isr_day_night.mat');
        
        
    end;
    if isunix == 1,
        load /data/backup/mnair/ace_tensor/acedata/ace_2000_2010.mat ace_all;;
        ace_all(:,1) = ace_all(:,1)+(17)/(60*24);%Try to advance 60 VS model/ ace time minutes
        %load '/data/backup/mnair/ace/jcamarca_isr_day_night.mat';
        load /data/backup/mnair/ace/jcamarca_isr_fejer.mat eef;
% eef data contains        
% eef(:,1) = fday - datenum(2000,1,1);
% eef(:,4) = fday_lt;
% eef(:,6) = drift*23000/1e6;
% eef(:,15) = fejer(:,2)*23000/1e6;
%readme = 'eef data contains r  eef(:,1) = fday - datenum(2000,1,1);  eef(:,4) = fday_lt; ... eef(:,6) = drift*23000/1e6; ... eef(:,15) = fejer(:,2)*23000/1e6'
 
        % To use the existing codes, just fill the data to the variables
        fday = [ eef(:,1) + datenum(2000,1,1)]';
        drift = [ eef(:,6) - eef(:,15) ]';
    end;
    
end;

options.Eeff_cut_off = 8; % The cut off amplitude for IEF .LUHR & MAUS EPS 2010 used a value of 8
options.ace_interp_method = 'spline'; %FIXED spline
options.des_int = 0.25 ;
options.ap_lower_limit = 20;
load /data/backup/mnair/longp/aplist.mat;


% Find data for each hour in a day
% the idea is to find cross correlation coefficient between the
% IEF Ey segments and the EEF data pieces centered around each LT hour.
% Hope to see the +ve day and -ve night correlation
% Use the cc and some division / scalinf factor to give LT dependency to
% the TF


% remove large values in ACE data

L = ace_all(:,2) < 0;
temp = abs(ace_all(:,2));
Eeff = options.Eeff_cut_off * temp ./ sqrt (options.Eeff_cut_off^2 + temp.^2);
ace_all(:,2) = Eeff;
ace_all(L,2) = ace_all(L,2) * -1;

% interpolate and down sample

L = isnan(ace_all(:,2));
y = interp1(ace_all(~L,1),ace_all(~L,2),ace_all(:,1),options.ace_interp_method);
ace_down = resample(y,1,options.des_int*60/5); % Checked the resampling time axis OK
ace_inter = interp1(ace_all(1:options.des_int*60/5:end,1),ace_down,ace_all(:,1), options.ace_interp_method);
ace_all(:,2) = ace_inter;



%%

n_ut = 1;
ut_interval = 0.5;%hours
window_length = 3; % hours
min_piece_length = ( window_length * 60 / 5 ) -  1;
xplot = [];

datacount = zeros([1, (24 / ut_interval) + 1 ]);
correlation_ut = zeros([1, (24 / ut_interval) + 1]);
correlation_sig = zeros([1, (24 / ut_interval) + 1]);
correlation_rlo = zeros([1, (24 / ut_interval) + 1]);
correlation_rup = zeros([1, (24 / ut_interval) + 1]);
data_ratio = zeros([1, (24 / ut_interval) + 1]);
ratio_count = zeros([1, (24 / ut_interval) + 1]);
corr_error = data_ratio;
for i = 0 : ut_interval: ( 24 ),
    
    
    
    if ( i + window_length) >= 24,
        K = [];
        A = [];
        upper_limit = i + window_length - 24;
        L = (fday - floor(fday) >= (i/24)) | (fday - floor(fday) < upper_limit /24 );
        this_ut_fday = fday(L);
        this_ut_drift = drift(L);
        A = find(diff(this_ut_fday) > 0.6);
        K(:,1) = [1 A+1];
        K(:,2) = [A length(this_ut_fday)];
        
        
    else
        L = (fday - floor(fday) >= (i/24)) & (fday - floor(fday) < ( (i + window_length) /24));
        
        this_ut_fday = fday(L);
        this_ut_drift = drift(L);
        K = floor(this_ut_fday);
        A = unique(K);
        
    end;
    
    
    
    
    
    
    for j = 1 : length(A),
        
        if ( i + window_length) >= 24,
            
        
                piece_of_drift = this_ut_drift(K(j,1):K(j,2));
                piece_of_fday =  this_ut_fday(K(j,1):K(j,2));
        
        else
            
            L = floor(this_ut_fday) == A(j);
            piece_of_drift =  this_ut_drift(L);
            piece_of_fday =  this_ut_fday(L);
        end;
        
        
        if length(piece_of_fday) > 10,
            
            
            
            kkk = piece_of_fday(1):5/1440:piece_of_fday(end);
            
            %fprintf('%d %d %d\n', i, length(this_ut_drift),length(piece_of_fday));
            
            %get ace data for this piece of data
            
            L = ace_all(:,1) >= ( piece_of_fday(1) - 10/(60*24) ) & ...
                ace_all(:,1) <= ( piece_of_fday(end) + 10/(60*24) );
            
            ace_time = ace_all(L,1);
            ace_data = ace_all(L,2);
            
            L = fday_ap >=  floor(piece_of_fday(1)) & fday_ap <=  ceil(piece_of_fday(end));
            mean_ap = mean(ap(L));
            
            piece_of_ace = interp1(ace_time, ace_data, kkk);
            piece_of_drift = interp1(piece_of_fday, piece_of_drift, kkk);
            
            %fprintf('%d %d %d %d %f \n', i, length(this_ut_drift),length(piece_of_fday), length (piece_of_ace), mean_ap );
            
            if ~any(isnan(piece_of_ace))   && ~any(isnan(piece_of_drift)) && ...
                    length(piece_of_ace)   >=  min_piece_length  && ...
                    length(piece_of_drift) >=  min_piece_length && ...
                    mean_ap > options.ap_lower_limit,
                
                %[R,P,RLO,RUP]=corrcoef(...) also returns matrices RLO and RUP,
                %of the same size as R, containing lower and upper bounds for a 95%
                %confidence interval for each coefficient.
                
                [a,p,rlo,rup] = corrcoef(piece_of_ace, piece_of_drift);
                datacount(n_ut) = datacount(n_ut )  + 1;
                correlation_ut(n_ut) = correlation_ut(n_ut) + a(2,1);
                correlation_mat(n_ut,datacount(n_ut)) = a(2,1);
                correlation_sig(n_ut) = correlation_sig(n_ut) + p(2,1);
                correlation_rlo(n_ut) = correlation_rlo(n_ut) + rlo(2,1);
                correlation_rup(n_ut) = correlation_rup(n_ut) + rup(2,1);
                corr_error(n_ut) = corr_error(n_ut) + (1 - (a(2,1))^2) / sqrt( length(piece_of_ace ) - 1 );
                
                if ~any(piece_of_ace == 0),
                    %data_ratio(n_ut) = data_ratio(n_ut) + sum((piece_of_drift./piece_of_ace).^2);
                    data_ratio(n_ut) = data_ratio(n_ut) + sum(abs(piece_of_drift))/sum(abs(piece_of_ace));
                    ratio_count(n_ut) = ratio_count(n_ut) + 1;
                end;
                
                
            end;
            
        end;
        
    end;
    xplot(n_ut) = i + window_length  /2;
    n_ut = n_ut + 1;
    
end;


%% Alternative way to calculate correlation
% THERE is an un identified issue with this script. Use the above one
% Spent almost a day on this

%%
% n_ut = 1;
% ut_interval = 0.5;%hours
% window_length = 2; % hours
% min_piece_length = ( window_length * 60 / 5 ) -  1;
% 
% datacount = zeros([1, (24 / ut_interval) + 1 ]);
% correlation_ut = zeros([1, (24 / ut_interval) + 1]);
% correlation_sig = zeros([1, (24 / ut_interval) + 1]);
% correlation_rlo = zeros([1, (24 / ut_interval) + 1]);
% correlation_rup = zeros([1, (24 / ut_interval) + 1]);
% data_ratio = zeros([1, (24 / ut_interval) + 1]);
% ratio_count = zeros([1, (24 / ut_interval) + 1]);
% corr_error = data_ratio;
% %%
% 
% data_index = [];
% n_ut = 1;
% for j = 0 : ut_interval: 24,
%     %for j = 0.5 :0.5,
%     
%     
%     nd = 1;
%     np = 0;
%     
%     
%     
%     for i = 1: length(fday) - 1,
%         
%         if j + window_length >= 24 & (fday(i) - floor(fday(i)) )*24 <= window_length,
%             upper_limit = (fday(i) - floor(fday(i)) )*24 + 24;
%             %         display(i);
%             %         pause;
%         else,
%             upper_limit = (fday(i) - floor(fday(i)) )*24;
%         end;
%         
%         
%         if (fday(i+1) - fday(i) <= 6/1440 & upper_limit >= j...
%                 & upper_limit <= j + window_length + 5/60 )
%             
%             
%             np = np + 1;
%             
%             
%             
%             if (i - np ) > 0 && abs( fday(i) - fday(i - np) - window_length/24 ) < 1/1440,
%                 % display(np);
%                 
%                 data_index(n_ut,nd,2) = i;
%                 data_index(n_ut,nd,1) = i - np;
%                 np = 0;
%                 nd = nd + 1;
%                 
%             end;
%             
%         else
%             np=0;
%             
%             
%         end;
%         
%     end;
%     
%     
%     
%     n_ut = n_ut + 1;
% end;
% 
% n_ut = 1;
% for i = 1:size(data_index,1),
%     
%     for j = 1:size(data_index,2),
%         
%         
%         
%         if      data_index(i,j,2) - data_index(i,j,1) > 1,
%             piece_of_drift =  drift(data_index(i,j,1):data_index(i,j,2));
%             piece_of_fday =    fday(data_index(i,j,1):data_index(i,j,2));
%             
%             %fprintf('%d %d %d\n', i, length(this_ut_drift),length(piece_of_fday));
%             
%             %get ace data for this piece of data
%             
%             L = ace_all(:,1) >= ( piece_of_fday(1) - 10/(60*24) ) & ...
%                 ace_all(:,1) <= ( piece_of_fday(end) + 10/(60*24) );
%             
%             ace_time = ace_all(L,1);
%             ace_data = ace_all(L,2);
%             
%             L = fday_ap >=  floor(piece_of_fday(1)) & fday_ap <=  ceil(piece_of_fday(end));
%             mean_ap = mean(ap(L));
%             
%             % mean_ap = 40;
%             
%             piece_of_ace = interp1(ace_time, ace_data, piece_of_fday);
%             
%             %fprintf('%d %d %d %d %f \n', i, length(this_ut_drift),length(piece_of_fday), length (piece_of_ace), mean_ap );
%             
%             if ~any(isnan(piece_of_ace))   && ~any(isnan(piece_of_drift)) && ...
%                     length(piece_of_ace)   >=  min_piece_length  && ...
%                     length(piece_of_drift) >=  min_piece_length && ...
%                     mean_ap > options.ap_lower_limit,
%                 
%                 %[R,P,RLO,RUP]=corrcoef(...) also returns matrices RLO and RUP,
%                 %of the same size as R, containing lower and upper bounds for a 95%
%                 %confidence interval for each coefficient.
%                 
%                 [a,p,rlo,rup] = corrcoef(piece_of_ace, piece_of_drift);
%                 datacount(n_ut) = datacount(n_ut )  + 1;
%                 correlation_ut(n_ut) = correlation_ut(n_ut) + a(2,1);
%                 correlation_mat(n_ut,datacount(n_ut)) = a(2,1);
%                 correlation_sig(n_ut) = correlation_sig(n_ut) + p(2,1);
%                 correlation_rlo(n_ut) = correlation_rlo(n_ut) + rlo(2,1);
%                 correlation_rup(n_ut) = correlation_rup(n_ut) + rup(2,1);
%                 corr_error(n_ut) = corr_error(n_ut) + (1 - (a(2,1))^2) / sqrt( length(piece_of_ace ) - 1 );
%                 
%                 if ~any(piece_of_ace == 0),
%                     %data_ratio(n_ut) = data_ratio(n_ut) + sum((piece_of_drift./piece_of_ace).^2);
%                     data_ratio(n_ut) = data_ratio(n_ut) + sum(abs(piece_of_drift))/sum(abs(piece_of_ace));
%                     ratio_count(n_ut) = ratio_count(n_ut) + 1;
%                 end;
%                 
%                 
%             end;
%         end;
%         
%     end;
%     
%     xplot(n_ut) = i + window_length  /2;
%     n_ut = n_ut + 1;
%     
% end;
%% plot the correlation
%xplot = ut_interval/2 : ut_interval : 24;

% find the median correlation

for i = 1:49,
L = correlation_mat(i,:) ~= 0;
corrmed_ut(i) = median( correlation_mat(i,L));
corr_rb(i) =  robustfit(zeros(sum(L),0),correlation_mat(i,L));
end;

% find the standarad deviation
for i = 1:49,
L = correlation_mat(i,:) ~= 0;
sstd(i) = std(correlation_mat(i,L));
end;

%plot all the data

for i = 1:49,
L = correlation_mat(i,:) ~= 0;
dummy =  correlation_mat(i,L) ;
for j = 1:length(dummy),
plot(xplot(i),dummy(j),'r*');
hold on;
end;
end;

%Generate the time axis
xplot = (0 : ut_interval : 24) + window_length/2;

% Which are GT 24
L = xplot >= 24;

% circile shift the data > 24 as the first
xplot = circshift(xplot,[1,sum(L)]);
corplot = circshift(correlation_ut,[1,sum(L)]);
dataplot = circshift(datacount,[1,sum(L)]);
corerr = circshift(corr_error,[1,sum(L)]);
corsig = circshift(correlation_sig,[1,sum(L)]);
corerrlow = circshift(correlation_rlo,[1,sum(L)]);
corerrup = circshift(correlation_rup,[1,sum(L)]);
corrmed = circshift(corrmed_ut,[1,sum(L)]);
sstd = circshift(sstd,[1,sum(L)]);
               
% Reduce 24 from LT > 24
L = xplot >= 24;

xplot(L) = xplot(L) - 24;
 
% Any duplicate entries ?
[y,ia] = unique(xplot);
  
if any(ia)
     
xplot = xplot(ia);
corplot = corplot(ia);
dataplot = dataplot(ia);
corerr = corerr(ia);
corsig = corsig(ia);
sstd = sstd(ia);

corerrlow = corerrlow(ia);
corerrup = corerrup(ia);
corrmed = corrmed(ia);
end;
 

%% plot the data
figure(3);
 
%errorbar(xplot, corplot./dataplot,corerr./dataplot, corerr./dataplot ,'ks');

% Using P test as the correlation significance
%errorbar(xplot, corplot./dataplot, (corsig./dataplot)/2, (corsig./dataplot)/2 ,'ks');
%errorbar(xplot, corplot./dataplot, (corerrlow./dataplot), (corerrup./dataplot) ,'ks');
errorbar(xplot, corplot./dataplot, sstd/2, sstd/2 ,'ks');
%errorbar(xplot, corrmed, (corerrup./dataplot), (corerrlow./dataplot) ,'ks');

xticks = 0:2:24 ;

set(gca,'XTick', xticks);


xticks_lt = xticks - 5;

xticks_lt (xticks_lt <= 0) = xticks_lt (xticks_lt <= 0) + 24;

set(gca,'XTickLabel',  reshape(sprintf('%2.0f',xticks_lt),[2,length(xticks_lt)])' );

set(gca,'FontSize',16);
xlabel('Local Time (Hours)');
ylabel('Correlation');
legend('Jicamarca ISR - ACE IEF Ey');
grid on;
hold on

% plot the polyfit

% xplot_polyfit = [-1.5 -1.0 -0.5 xplot 24 24.5 25];
% 
% [p,s ] = polyfit(xplot_polyfit ,...
%     [corplot(end-2:1:end) corplot corplot(1:3) ]./ ...
%     [dataplot(end-2:1:end) dataplot dataplot(1:3) ], 7);

xplot_polyfit = [-fliplr(xplot)-0.5 xplot xplot+24];

%plot([-fliplr(xplot)-0.5 xplot xplot+24] ,[corplot corplot corplot ]./[dataplot dataplot dataplot ])

[p,s ] = polyfit(xplot_polyfit ,[corplot corplot corplot ] ./ [dataplot dataplot dataplot ], 21);

 

% plot(xplot_polyfit ,...
%     [corplot(end-2:1:end) corplot corplot(1:3) ]./ ...
%     [dataplot(end-2:1:end) dataplot dataplot(1:3) ]);


%
y = polyval(p,xplot);
plot(xplot, y,'r','LineWidth',2);
axis([-1,25,-.8,.8]);

%% make a polyfit on LT axis for pp_lt_response.m program
%
%factor =  2.6537 ; %1/max(y);

factor = 1/max(y);

correlation_data = [corplot./dataplot corplot./dataplot corplot./dataplot] *  factor ;

%  factor (=1./mean(correlation(14:21)) for scaling the data.
%  the idea is the the mean power between 9-16LT should have a
%  multiplication factor 1.

local_t_xscale = [-24.0 : 0.5 : 47.5] - 5;

%  [s] = csaps(local_t_xscale,correlation_data);
%  sn = fnxtr(s);%This methods gives superior extrapolation !
%  correlation_smoothed = ppual(sn,local_t_xscale);
%  correlation_smoothed = correlation_smoothed(25:48);

%  correlation_smoothed = v(25:48);


[p,s ] = polyfit(local_t_xscale ,correlation_data, 21);

y = polyval(p,0:1:24);
plot(0:1:24, y,'r','LineWidth',2);
axis([-1,25,-1.1,1.1]);


% may 3, 2012. While above LT scaling has been used in the paper, the
% implementation model uses slightly different scaling (two reasons : all
% data and an inferior polyfit). This need to be fixed.
